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1.  Introduction 

Solutions  of  advection-dispersion  equation  (ADE)  may  be 
used  to  predict  the  concentration  of  solutes  in  unsteady 
groundwater  flow.  Advection  causes  the  contaminant  plum 
to  flow  in  the  direction  of  groundwater  water  without  any 
change  in  the  shape.  Dispersion  of  plum  arises  due  to  the 
variation  in  groundwater  velocity.  The  heterogeneity  of  the 
porous  medium  is  responsible  for  dispersion.  The  solute 
transport  in  heterogeneous  aquifer  is  thus  the  combined 
process  of  advection  and  dispersion.  The  dispersion  in 
the  direction  of  groundwater  flow  is  called  longitudinal 
dispersion  and  the  transverse  dispersion  is  perpendicular  to 
the  groundwater  flow  direction.  The  ADE  can  be  derived 
using  Fick’s  law  and  low  of  conservation  of  mass. 

Analytical  solutions  in  one-dimensional  problems 
through  semi-infinite  or  finite  porous  media  have  been 
presented  by  several  researchers:  (Mazaheri  et.  al.  2013, 
Kumar  et.  al.  2010,  Marino  et  al.  19 74,  Singh  et  al.  2008) 
etc.  The  objective  of  this  work  is  to  derive  an  approximate 
analytical  solution  of  ADE  with  the  help  of  Optimal 
Homotopy  Analysis  Method  (OHAM).  In  this  work,  an 
approximate  analytical  solution  of  one-dimensional  ADE  in 
heterogeneous  finite  aquifer  is  derived  for  continuous  time 
dependent  input  source  concentration  of  increasing  nature. 


2.  Problem  Formulation 

We  assume  that  the  solute  transport  is  primarily  one¬ 
dimensional.  The  solutes  are  assumed  to  be  conservative 
in  the  unsteady  field.  Consider  an  isotopic  heterogeneous 
finite  aquifer  of  length  L=1  km. 

Let  C{x,t)  be  the  contaminant  concentration  in  the 
aquifer,  u  the  groundwater  velocity  and  DL  the  longitudinal 
dispersion  coefficient.  Then  the  partial  differential  equation 
describing  the  one  dimensional  advection  and  dispersion 
is  [ii] 


dC  .  d  ,  v  d 
(uC)  = 


dt  dx 


dx 


A 


dC 

dx 


(1) 


Expanding  (1)  we  obtain 

3Dl 


dC 

dt 


dx  , 


d£<C—  =  D  — 

dx  dx  L  dx2 


(2) 


In  heterogeneous  medium,  porosity  changes  with  position. 
The  velocity  is  non  uniform  as  it  depends  on  porosity.  Hence 
the  velocity  of  the  flow  field  transporting  the  contaminants 
is  considered  spatially  dependent.  Let  velocity  at  the  origin 
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x  =  0  of  the  domain  be  u( 0  which  increases  to  u0  (1  +  b)  at 
x  =  L  where  a  real  constant  b<l  ensures  that  the  change  in 
velocity  is  of  small  order  i.e.  the  laminar  condition  of  the 
flow  is  not  affected  (Kumar  et  al. ,  2010)  .  The  expression 
for  velocity  at  any  position  x  may  be  linearly  interpolated  as 

,  .  x  —  L  x  —  0 

u{x)=-^—^u0+j—^u0{\.  +  b)  (3) 


The  concentration  change  is  very  negligible  at  the  other 
end  L=  1  km  of  the  aquifer.  So  the  solute  transport  may 
not  be  affected  at  the  other  end  of  the  aquifer.  Hence  we 
prescribe  the  outlet  boundary  condition  as  no  flux  boundary 
condition,  i.e. 

dC 

—  (l,f)  =  0,f>0  (8) 

ox 


Simplifying  we  have 


u(x)  =  u0  (1  +  ax)  (4) 

b 

where  a  =  —  is  a  constant  less  than  1  and  serves  as  a 
L 

heterogeneity  parameter.  Its  different  values  represent  media 
of  varying  heterogeneity. 

Since  mechanical  dispersion  depends  on  the  flow, 
it  is  expected  to  increase  with  increasing  flow  speed.  The 
dispersion  parameter  is  considered  to  be  proportional  to  the 
square  of  the  velocity. 

Dl  =  Dq  (1  +  ax)2  (5) 

where  DQ  =  au2  where  a  is  known  as  longitudinal 
dispersivity. 

Using  equation  (4)  and  equation  (5)  in  equation  (2), 
we  obtain 


We  solve  equation  (6)  with  boundary  conditions  equation 
(7)  and  equation  (8)  using  optimal  homotopy  analysis 
method. 

3.  Solution  of  the  Problem  Using  OHAM 

To  solve  the  problem  by  OHAM,  we  choose  the  initial  guess 
of  the  solution  C(x,  t)  as 


C0(x,t)  =  t(e  x  -\- xe  1 ) 


(9) 


which  satisfies  boundary  conditions  equation  (7)  and 
equation  (8).  The  linear  operator  is  selected  as 


L  [</>(x,t;q)\ 


d2<j)(x,t;q) 

d? 


(10) 


satisfying  the  property 


L  [f]  =  0  when  f  =  0. 


(ID 


dC  ,  „  ^  ,dC 

— — F  («0  —  2aD0)(l  +  ax)—— 
ot  ox 

d2C 


-\~aUqC  —  D0  (1 H-  ax) 


dx 2 


(6)  where  0  <  q  <  1  is  the  embedding  parameter. 

Further  according  to  equation  (6),  the  other  operator 
is  defined  as 


which  is  linear  advection-dispersion  equation  with  variable 
coefficients.  The  solution  C{x,t)  of  this  equation  represents 
the  concentration  of  contaminants  at  any  position  and  at 
any  time. 

The  boundary  condition  at  the  origin  is  of  uniform 
nature  as  defined  by  many  researchers.  It  means  that  the  input 
concentration  at  the  origin  of  the  medium  (water  bodies 
on  the  surface,  groundwater  level)  and  hence  its  source  of 
pollution  on  the  earth  surface  remains  uniform  at  all  times. 
But  this  is  not  the  real  situation  at  all  times.  Actually,  the 
pollution  at  the  source  and  so  the  input  will  increase  with 
time  due  to  increasing  human  activities  on  the  surface.  So  we 
have  considered  the  inlet  boundary  condition  as 

C(0,t)  =  t  ,  t>0  (7) 


N[0(#,  t;  q)\  =  D0  (1  +  ax)2  — ^  —  —  au0(j)(x,  t;  q) 

dx 

-H-2  aD^l+ax)aJ^l  (12) 

dx 

d(j)(x,t',q) 

dt 

Let  cQ  denote  a  nonzero  auxiliary  parameter.  According  to 
Liao  [iii] ,  the  zeroth  order  deformation  equation  is 

(1  -  q) L[<f)(x,  t\  q)  -  C0  (x,  t)\  =  c0qH(x ,  t) N[0(x,  t;  q)}  (13) 

where  H(x,t)  is  non-zero  auxiliary  function  and  <f>(x,t;q) 
is  an  unknown  function. 


ISSN  No.:  2278-9561  (Print)  ISSN  No.:  2278-957X  (Online)  Registration  No. :  CHAENG/20 13/49583 


Math.  J.  Interdiscip.  Sci..Vol.  7,  No.l,  Sep.  2018 


pp.17 


Clearly,  when  q  =  0  and  q  =  1  ,  we  have  from  equation 
(11)  and  equation  (13), 


0)  =  C0(x,t) 


(14) 


Differentiating  equation  (13)  m  times  with  respect  to  q  and 
then  dividing  them  by  ml  and  then  taking  q  =  0  ,  we  obtain 
high  order  deformation  equations 


L  [Cm{x,t)-XmCm-^x^)\  =  cQH(x,t)Rm(x,t)  (19) 


and  (j)(x,t;l)  =  C(x,t) 


(15) 


subject  to  the  boundary  conditions 


Hence  as  q  increases  from  0  to  7,  the  solution  c/)(x,t;q ) 
deforms  from  the  initial  choice  C0  (x,  t)  to  the  exact  solution 
C(x,t)  of  equation  (6).  The  determination  of  <fi(x,t;q ) 
depends  on  the  proper  choices  of  the  operator  Z,  the  initial 
choice  C0  (x,  t)  and  the  parameter  c0 .  We  assume  that  all  of 
them  are  properly  chosen,  the  Maclaurin  series 


oo 

4>{x,t\q)  =  C0{x,t)  +  'Y^Cm{x,t)qm  (16) 

m= 1 

exists  and  converges  at  q  =  1  .  Hence  we  have  the  homotopy 
series  solution 


c{x,t)=c0(x,t)+ y yjm  (x,  t )  (17) 

m= 1 


where  Cm(x,t)  = 


1  dm<p(x,t,q ) 

~m\  frf 


l?-o 


(18) 


C(0,t)  =  0, 


dCm 

dx 


(l,f)  =  0,  m>  1 


where 


RJx’t)  =  -, - — - - I, 


dqm~l  '?=0 
d2C 

u  V  i 


—  ax) 


dx 1 


— (u0  —  2aD0 )  (1  +  ax) 

dCm Y 


dCm_, 

dx 


aUvPm-\ 


dt 


and  Xm  ~ 


0  if  m  <  1, 
1  if  m  >  1 . 


(20) 


(21) 


(22) 


We  assume  H(x,t)  =  1  for  simplicity.  Equations  (19)  are 
nonhomogeneous  linear  ordinary  differential  equations 
with  constant  coefficients  for  all  m  >  1  .  Solving  equations 
(19)-(20)  for  m  =  1 ,  we  obtain  the  first  order  homotopy 
approximation  as 


Cl(x,t)  =  c0 


cl  D^tx  e  x  H-  (2a  D0  T  2 T  clu 0  j  *  -\-(uq  -f  2^  D0  T  2 T  D0  T  clu 0  j  ^ 


3  -1 

x  e 


-  + 


aD{)e 


tx2  + 


a2  D0e  1 


a 

/x - 


-  +  1 


-f  ( Dr,e  2ur>e  T  2aUr]e  )tx  —  ( 2 a  T  2uDr)  -f  Dn  T  T  )  t 


(23) 


Similarly  the  mth  order  homotopy  approximation  C  (x,t) 
can  be  obtained  for  successive  values  of  m.  The  convergence 
of  the  mth  order  homotopy  approximation  depends  on  the 
proper  choice  of  c0 .  The  homotopy  series  solution  can  then 
be  written  as 


residual  at  the  mth  order  homotopy  approximation  denoted 
by  Em  and  defined  by  (Liao  et  al.  20 1 0)  as 


- £ - EE 

(M  +  l)(7V  +  l)t^ 


N 


Yfm 


m’n 


(25) 


C(x,t)  =  C0(x,t)  +  Cl(x,t)  + ...  (24) 

The  solution  represents  solute  or  contaminant  concentration 
at  distance  x  for  any  time  t  whose  convergence  depends  on 
the  proper  choice  of  the  convergence-control  parameter  c0. 
To  select  proper  value  of  c0.,  we  use  the  discrete  squared 


As  fast  the  residual  Em  decreases,  the  accuracy  of  the 
corresponding  homotopy  approximation  increases.  At  the 
mth  order  approximation,  the  corresponding  optimal  value 
of  the  convergence-control  parameter  c0  is  given  by  the 
minimum  of  Em  corresponding  to  a  nonlinear  algebraic 
equation  to  be  solved  from 
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dEm(c0) 

dCr, 


=  0 


(26) 


This  approach  for  obtaining  the  optimal  value  of  c0  has  been 
used  for  solving  a  number  of  problems  for  nonlinear  ordinary 
and  partial  differential  equations  by  many  researchers  (Liao, 

SquaredResidual 


2010,  Prajapati  et  al.  2015,  2016,  2017,  Liao  et  al.  2012, 
2013,  vejravelu  et  al.,  2012).  We  have  obtained  the  solution 
up  to  10th  order  approximation.  The  optimal  value  of  CQ 
is  found  by  the  minimum  of  El0  using  Mathematica. 
Here  E10  attains  its  minimum  value  5.23828  xlCT7  at 
c0  =  —2.42025  which  we  can  observe  in  Figure  1  also.  We 
take  M  =  N  =  50  for  finding  Em  . 


Figl  The  discrete  squared  residual  at  10th  order  homotopy  approximation 


Table  1.  The  discrete  squared  residual  of  governing  equation  (6)  by 
means  of  c0  =  —2.42025  • 


Order  of  approximation  m 

Discrete  Squared  Residual  E 

2 

1.24554E-1 

4 

3.53453E-3 

6 

2.38709E-4 

8 

9.68075E-6 

10 

4.7776E-7 

From  the  Table  1,  we  can  see  that  the  squared  residual 
decreases  as  we  increase  the  order  of  approximation.  So  the 
10th  order  homotopy  approximation  is  accurate. 

4.  Numerical  and  Graphical  Representation 

BVPh,  a  Mathematica  package,  is  used  to  obtain  numerical 
values  of  concentration.  Table  2  indicates  the  numerical 
values  of  the  solute  concentration 


C(x,t) 


kg 


[kmf 


for  different  distance  x  (km)  and  time  t  (year)  up  to  10th 
order  approximation  using  c0  =  —2.42025  . 

We  have  considered 


L  =  \  (km) ,  a  =  Q.\(km)  ' ,  u0  =0.11 
(km) 


E>0  =0.21 


yr 


km 


yr 


and 


to  obtain  numerical  values  of  the  solution. 

Figure  2  shows  the  concentration  profiles  of  the 
contaminant  for  fixed  values  of  times  .  Numerical  values  of 
Table  2  are  used  for  Figure  2.  Figure  3  shows  the  graph  of 
concentration  for  fixed  positions.  We  use  numerical  values 
of  Table  2  for  Figure  3  also. 


Table  2.  Numerical  values  of  the  solute  concentration  C(x,  t). 


X 

t=l 

t=l.l 

t  =1.2 

t  =1.3 

t  =1.4 

t=1.5 

t  =1.6 

t  =1.7 

0.0 

1.00000 

1.10000 

1.20000 

1.30000 

1.40000 

1.50000 

1.60000 

1.70000 

0.01 

0.96338 

1.06334 

1.16330 

1.26326 

1.36322 

1.46318 

1.56314 

1.66310 

0.02 

0.92713 

1.02705 

1.12697 

1.22689 

1.32681 

1.42673 

1.52665 

1.62657 

0.03 

0.89124 

0.99112 

1.09100 

1.19088 

1.29076 

1.39064 

1.49052 

1.59040 
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0.04 

0.85571 

0.95555 

1.05539 

1.15523 

1.25507 

1.35491 

1.45475 

1.55459 

0.05 

0.82054 

0.92034 

1.02014 

1.11995 

1.21975 

1.31955 

1.41935 

1.51915 

0.06 

0.78574 

0.88550 

0.98526 

1.08503 

1.18479 

1.28455 

1.38432 

1.48408 

0.07 

0.75129 

0.85102 

0.95074 

1.05047 

1.15019 

1.24992 

1.34965 

1.44937 

0.08 

0.71721 

0.81690 

0.91659 

1.01628 

1.11596 

1.21565 

1.31534 

1.41503 

0.09 

0.68349 

0.78314 

0.88280 

0.98244 

1.08210 

1.18175 

1.28140 

1.38105 

0.10 

0.65014 

0.74975 

0.84936 

0.94898 

1.04859 

1.14821 

1.24782 

1.34743 

Ctr,  t) 


Fig  2  The  spatial  distribution  of  solute  concentration  for  fixed  values  of  times 
C(x,  t ) 


Fig  3  The  temporal  distribution  of  solute  concentration  for  fixed  values  of  distances 


5.  Conclusion 

The  approximate  analytical  solution  of  advection-dispersion 
equation  with  variable  coefficients  is  obtained  by  optimal 
homotopy  analysis  method  with  time-dependent  input 
source  concentration.  The  contaminant  transport  behaves  as 


expected  i.e.  the  concentration  at  a  fixed  time  decreases  as  the 
distance  increases  and  the  concentration  at  a  fixed  position 
advances  with  time.  The  solution  is  useful  as  a  preliminary 
predictive  tool  for  simulating  the  solute  migration  in  aquifer 
due  to  the  release  of  a  time-dependent  source. 
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